OSM Beaches

We are going to use a beach shapefile from OpenStreetMaps and calculate the length and the width of the identified beaches within the maps. We will also need to link these to Hawaii poylgons like what was done in the buffer approach. Here what we do is use a buffer to link houses that adjoin the polygon beaches

Code

beaches <- read_sf("/Users/ashleylowe/Trails/beaches_meters.shp")
beaches <- st_transform(beaches, 32605)
beaches_buffered <- st_buffer(beaches, dist = 10)
beaches_buffered <- st_transform(beaches_buffered, 4326)
source("_common.R")



intersections <- st_intersects(HI_2023_shapefile, beaches_buffered)
HI_intersections <- HI_2023_shapefile[lengths(intersections) > 0, ]


ggplot() +
    geom_sf(data = HI_intersections, fill = "brown", ) +
    geom_sf(data = beaches_buffered, fill = "lightblue") +
  geom_sf(data =coastline, color = "lightgrey") +
    coord_sf(
        xlim = c(-158.2, -157.6),  # longitudes
        ylim = c(21.2, 21.6),        # latitudes in correct order
        expand = FALSE
    ) +
    theme_minimal() +
    labs(title = "Buffered Beaches (Blue) and Adjoining Properties (Brown)")

Using the Adjoining properties of these beaches we will now merge the beach names in OSRM (note not all beaches have a name), beach width and length based on the polygon calculate using the skeleton of the polygon in QGIS.

Code

# Create a vector of indices for the first overlapping shp2 feature per shp1 row
HI_2023_shapefile$beachname <- sapply(intersections, function(idx) {
  if (length(idx) == 0) {
    return(NA_character_)
  } else {
    paste(beaches_buffered$name[idx], collapse = ", ")
  }
})

# Create a vector of indices for the first overlapping shp2 feature per shp1 row
HI_2023_shapefile$BeachWidth <- sapply(intersections, function(idx) {
  if (length(idx) == 0) {
    return(NA_character_)
  } else {
    paste(beaches_buffered$BeachWidth[idx], collapse = ", ")
  }
})

HI_2023_shapefile$Beachlength <- sapply(intersections, function(idx) {
  if (length(idx) == 0) {
    return(NA_character_)
  } else {
    paste(beaches_buffered$LENGTH[idx], collapse = ", ")
  }
})

Coral Reefs

I am also interested in the coral reefs and if those are incorporated into housing prices. Using the [Asner et al](https://www.pnas.org/doi/10.1073/pnas.2017628117) Using the property file we isolate the potential coastal properties following

• Selecting the coastline shapefile

• Create a buffer of 10meters for coastline

• Select Hawaii properties intersecting the coastline

• Buffer properties 300 meters to grab the areas coral cover average.

The following file shows properties which returned a value using this method.

Code

asner <- read_csv("/Users/ashleylowe/Hawaii - Recreation/Island_Condition/Coral/asner_coastline50meters.csv")

asner$costline_property=1

asner=HI_2023_shapefile %>%
    semi_join(asner, by = "objectid") %>%       # keep only matched rows
    left_join(asner %>% select(objectid, coral_asner), by = "objectid")  # bring in coral_cover
# Transform to WGS84
asner <- st_transform(asner, 4326)

tmap_mode("view")
center_lon <- -157.45
center_lat <- 20.75
zoom_level <- 7

tm_shape(asner) +
  tm_polygons(
    "coral_asner",
    palette = "viridis",
    title = "Coral Reef Cover",
    popup.vars = c("Coral Reef Cover" = "coral_asner")
  ) +
  tm_basemap() +
  tm_view(set.view = c(center_lon, center_lat, zoom_level))

Matching these to the housing locations across the islands

Code

# Make sure both have the same CRS
st_crs(PB_HawaiiInfo_joined) == st_crs(asner)  

# If FALSE, transform one:
#PB_HawaiiInfo_joined <- st_transform(PB_HawaiiInfo_joined, st_crs(Special_Management_Areas))

# Spatial join: attach polygon attributes to points
PB_HawaiiInfo_joined<- st_join(PB_HawaiiInfo_joined, asner, join = st_within)